rm(list=ls(all=TRUE))

library(MASS)
library(mc2d)

n.MC <- 2000
pars.MC <- matrix(NA,nrow=n.MC,ncol=3)
mu.hat <- rep(NA,length.out=n.MC) 
delta  <- rep(NA,length.out=n.MC)

n.sim <- 5000  ## sample size for each Monte Carlo simulation ##

X.mu  <- c(.5, 1.5)
X.sig <- rbind(c(.1,.1),c(.1,.2))

lamT1.act <- 0.7
lamL1.act <- 0.2
lamT0.act <- 0.9
Beta.act  <- c(-1,2,-1.3)

for (k in 1:n.MC) {

X<-mvrnorm(n.sim,X.mu,X.sig)
X<-cbind(1,X)

mu.act <- (1+exp(-X%*%Beta.act))^-1
theta  <- rbinom(n=n.sim,size=1,mu.act)

gam1 <- lamT0.act*(1-theta) + lamL1.act*theta
gam2 <- lamT1.act*theta
gam3 <- (1-lamT1.act-lamL1.act)*theta + (1-lamT0.act)*(1-theta)

Mprob <- cbind(gam1,gam2,gam3)

Y <- rmultinomial(n=n.sim,size=1,Mprob)
Y <- 1*Y[,1] + 2*Y[,2] + 3*Y[,3]

##remove observations from non-responders##

X <- X[Y<3,]
Y <- Y[Y<3]


## Make Y binary: 1 or 0 ##

Y <- 1*as.numeric(Y==2) + 0*as.numeric(Y==1)


## Set up function to calculate parameter values based on Newton-Raphson algorithm ##

NewRap <- function(par) {
	
Beta  <-par[1:length(par)]
mu    <- (1+exp(-X%*%Beta))^-1 

V <- mu*(1-mu)
X.til<-matrix(NA,nrow=nrow(X),ncol=ncol(X))
for (i in 1:length(Y)) {X.til[i,]<-V[i]*X[i,]} 

D <- Y - mu

Beta.j <- Beta + solve(t(X)%*%X.til,tol=.Machine$double.eps^2)%*%t(X)%*%D
Beta.j}

tol <- .Machine$double.eps^.5 # tolerance for stopping algorithm

M<-4000
K<-ncol(X)
init<-rep(0,K) # initial parameter estimates
pars.old <- NewRap(init)


for (j in 1:M){

pars.new <- NewRap(pars.old)

abs.diff <- abs(pars.new-pars.old)

if (max(abs.diff) < tol ) break

pars.old <- pars.new}

pars.MC[k,] <- pars.new

X2.l <- quantile(X[,3],.05)
X2.h <- quantile(X[,3],.95)

X.1 <- cbind(X[,1:2],X2.h)
X.0 <- cbind(X[,1:2],X2.l)

mu.ca   <- (1+exp(-X%*%pars.new))^-1
mu.1    <- (1+exp(-X.1%*%pars.new))^-1
mu.0    <- (1+exp(-X.0%*%pars.new))^-1 

mu.hat[k]<-mean(mu.ca)
delta[k]<-mean(mu.1-mu.0)

}

beta0.avg<-mean(pars.MC[,1])
beta1.avg<-mean(pars.MC[,2])
beta2.avg<-mean(pars.MC[,3])

beta0.bias <-beta0.avg-Beta.act[1]
beta1.bias <-beta1.avg-Beta.act[2]
beta2.bias <-beta2.avg-Beta.act[3]

beta0.mse <-mean((pars.MC[,1]-Beta.act[1])^2)
beta1.mse <-mean((pars.MC[,2]-Beta.act[2])^2)
beta2.mse <-mean((pars.MC[,3]-Beta.act[3])^2)

BETAS<-c(beta0.avg,beta1.avg,beta2.avg)
BIAS<-c(beta0.bias,beta1.bias,beta2.bias)
MSE<-c(beta0.mse,beta1.mse,beta2.mse)

QOIs<-cbind(mu.hat,delta)

save(QOIs,file="DirectQOIsLE5000.RData")

